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ABSTRACT 

Wc use recently published measurements of the kinematics, surface brightness and 
stellar mass-to-light ratio of the globular cluster NGC 2419 to examine the possibility 
that this Galactic halo satellite is embedded in a low-mass dark matter halo. NGC 2419 
is a promising target for such a study, since its extreme Galactocentric distance and 
large mass would have greatly facilitated the retention of dark matter. A Markov- Chain 
Monte Carlo approach is used to investigate composite dynamical models containing 
a stellar and a dark matter component. We find that it is unlikely that a significant 
amount of dark matter ( ^ 6% of the luminous mass inside the tidal limit of the 
cluster) can be present if the stars follow an anisotropic Michie model and the dark 
matter a double power law model. However, we find that more general models, derived 
using a new technique we have developed to compute non-parametric solutions to the 
spherical Jeans equation, suggest the presence of a significant dark matter fraction 
(approximately twice the stellar mass). Thus the presence of a dark matter halo around 
NGC 2419 cannot be fully ruled out at present, yet any dark matter within the 10' 
visible extent of the cluster must be highly concentrated and cannot exceed 1.1 x 
10^ Mq (99% confidence), in stark contrast to expectations for a plausible progenitor 
halo of this structure. 

Key words: globular clusters: individual (NGC 2419) — dark matter — stellar 
dynamics 



1 INTRODUCTION 

Dwarf spheroidal galaxies (dSphs) and globular clusters 
(GCs) have similar luminosities but are strikingly different 
in terms of their dark- matter (DM) content, as inferred from 
the kinematics of their stars: while dSphs (and their "ultra- 
faint" extension to low luminosity) galaxies are the most 
DM dominated systems, there is virtually no evidence that 
GCs have DM halos (Heggie & Hut 1996; Mashchenko & 
Sills 2005a; Bradford et al. 2011). This finding might repre- 
sent a problem for standard cold dark matter cosmology, in 
which GCs are expected to form within their own DM halos 
(Peebles 1984). However it is quite possible that GCs were 
originally embedded in massive DM halos that evolved dur- 
ing the cluster lifetime as a consequence of interaction with 
the cluster stars (Baumgardt & Mieske 2008) and stripping 
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by the tidal field of the host galaxy (Mashchenko & Sills 
2005b; Saitoh et al. 2006). Therefore, the central parts of 
GCs might be left relatively poor in DM, because the DM 
at the present time either populates the outer region of the 
cluster or has been stripped (Bekki & Yong 2012). However, 
this is unlikely to be the case in GCs observed to have stellar 
tidal tails, which argue against very extended massive DM 
components (Moore 1996). While it is not well established 
to what extent two-body relaxation and tidal stripping are 
effective in shaping hypothetical DM halos of GCs, it is clear 
that these effects are minimized in clusters with long two- 
body relaxation time and that orbit the outer parts of the 
host galaxy, where tidal effects are small. The prototype of 
these systems is the remote massive GC NGC 2419, located 
in the halo of the Milky Way at a distance of d = 87.5 kpc (Di 
Criscienzo et al. 2011) (Galactocentric distance of ~ 95 kpc), 
which has half-mass relaxation time ~ 43 Gyr (Harris 1996 
— 2010 edition). If NGC 2419 formed within a massive DM 
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halo, it is reasonable to expect that such a halo has main- 
tained its properties almost unaltered up to the present time. 

In this paper we use the recently published kinematic 
data set of 166 stars (Ibata et al. 2011a, hereafter Paper 
I) and the analysis of the stellar mass-to- light ratio Af,/L 
(Bellazzini et al. 2012) of the globular cluster NGC 2419 to 
study the DM content of this system. Paper I showed that 
the present kinematic and structural data on NGC 2419 can 
be fit acceptably well with models containing no DM; but 
this does not necessarily imply that models with DM cannot 
be as good. Previous studies have tried to put upper limits 
on the DM mass of NGC 2419 (Baumgardt et al. 2009; Con- 
roy, Loeb, & Spergel 2011). Baumgardt et al. (2009) have 
argued that the observed velocity dispersion profile of the 
cluster (based on 40 stars) is not consistent with a DM halo 
more massive than Mdm = lO^M© (assuming an isotropic 
stellar velocity distribution and NFW DM halo). Conroy, 
Loeb, & Spergel (2011) have put more stringent constraints 
(Mdm < IO^Mq), based on the observed outermost light 
profile of the cluster. However, the argument of Conroy, 
Loeb, & Spergel (2011) relies on the hypothesis that the 
cluster has relaxed via two-body relaxation, which is not 
thought to be the case for NGC 2419. 

Given the collisionless nature of NGC 2419 it is natural 
to explore, besides models with isotropic velocity distribu- 
tions, also models with radially or tangentially biased stel- 
lar orbits. Of course it is important to exclude very radially 
anisotropic systems, which are radial-orbit unstable. How- 
ever, in this respect it must be noted that stellar systems 
with massive DM halos can sustain relatively more radial 
kinetic energy than corresponding purely baryonic systems 
(Nipoti, Londrillo, & Ciotti 2002; Nipoti, Ciotti, & Londrillo 
2011). 

In Paper I we have shown that it is very problematic to 
explain the observed structure and kinematics of NGC 2419 
if modified Newtonian dynamics (MOND; Milgrom 1983) 
is the correct theory of gravity (but see Sanders 2012 and 
Ibata et al. 2011b). It must be stressed that this result does 
not necessarily imply that NGC 2419 does not contain a sub- 
stantial amount of DM. In fact, what we can infer on the ba- 
sis of the Ibata et al. (2011a) results is that it is unlikely that 
NGC 2419 has a MOND-equivalent DM halo, i.e. a DM halo 
with density distribution such that the overall gravitational 
potential of the cluster is the same as predicted by MOND 
(see Milgrom 1986; Nipoti, Ciotti, & Londrillo 2011, and ref- 
erences therein). Of course, there is no reason to expect that 
the putative DM halo of this globular cluster must mimic a 
MOND gravitational potential, so it is quite possible that 
there is a significant amount of DM in NGC 2419, but that 
this DM is not distributed as in a MOND-equivalent halo^. 

Throughout the paper we adopt the above distance to 
NGC 2419, which implies an angular scale of 25.452 pc per 
arcmin. Our conclusions are virtually independent of this 
specific choice, because the estimated uncertainty on the 
distance is small (~ 4%; Di Criscienzo et al. 2011). We 

^ It must be noted that stellar systems with MOND-equivalent 
halos can sustain more radial anisotropy in the central regions 
than MOND stellar systems (Nipoti, Ciotti, & Londrillo 2011), 
so, at least in principle, it is also possible that the cluster is well 
represented by very radially anisotropic models within a MOND- 
equivalent halo. 



note that our results on the dark-to-luminous mass ratio of 
NGC 2419, obtained for given d and stellar mass to light ra- 
tio (Mt/Lv), hold for all the combinations of d and M^/Lv 
such that (Mt/Lv) x d has the same value. In the present 
work we account for the uncertainty on M^/Lv, which is of 
the order of 30% (Bellazzini et al. 2012), so it is clear that 
our results are not affected by the smaller uncertainties on d. 
Other physical properties of the cluster are listed in Table 1 
of Paper I. 



2 MODELS WITHOUT DARK MATTER 

2.1 Distribution function based models 

Before analysing the case for the presence of dark matter in 
NGC 2419, it is instructive to consider models that possess 
only stars. In Paper I (see also Zocchi, Berlin, & Varri 2012), 
we presented several Michie model fits to the system in New- 
tonian gravity, and we showed that one of these models gives 
a statistically excellent representation of the cluster. How- 
ever, due to high computational cost, only a small grid of 11 
models were investigated in that contribution. To probe the 
parameter space much more finely, we re-wrote our Michie 
model code to allow it to be run on a parallel computer, 
thereby enabling us to calculate a large ensemble of Michie 
model fits, using a Markov-Chain Monte Carlo algorithm. 
This has the advantage of allowing us to properly incorpo- 
rate the measurement uncertainties and also yields posterior 
probability distribution functions for the fitted parameters. 
Furthermore, we are now able to use the stellar mass to light 
ratio measurement of Bellazzini et al. (2012) to constrain the 
solutions. 

The integration of the Michie models (optionally re- 
siding within within a dark matter potential) is outlined 
in Appendix A. We initially set the dark matter potential 
M^dm = so as to suppress that component, leaving 5 free 
fitting parameters: the Michie model stellar central density, 
the core radius, the anisotropy radius, the central potential 
and the stellar mass to light ratio^ . For a given set of these 5 
parameters, we construct a simulated surface brightness pro- 
file and the line of sight velocity distribution as a function of 
radius, which can be compared directly to the observations. 

As in Paper I, we convolve the model velocity distri- 
butions projected onto the line of sight by the estimated 
Gaussian errors associated to each kinematic datum i, to 
give a probability distribution fp{v, Ri) (we drop the model 
superscript j used in our earlier contribution, since we now 
consider a continuum of models defined by the above 5 pa- 
rameters). The likelihood of a model then becomes: 



\nL = J2^n[fp{v^,R^)]-J2 



25E2 



-+ln 



,(1) 



where the second summation over the (fc = 15) surface 
brightness profile data points is included to take account 
of the likelihood of the density profile, and the third term 



^ In Paper I, we analysed the effect on the line of sight velocity 
distribution of the presence of binaries in the cluster. The most 
likely binary fraction was found to be modest ( ^ 10%), so we 
decided to neglect this effect in the present contribution. 
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is the contribution from the observed stellar mass to light 
{Mt/Lv) ratio. 

In their analysis of recent WFC3 data of this cluster, 
Bellazzini et al. (2012) discuss in detail why the stellar mass 
to light ratio (Mt/Lv) must be constant, independent of 
radius. This is deduced from the uniformity of the luminos- 
ity function throughout the cluster, confirming an earlier 
analysis (Dalessandro et al. 2008) that reached the same 
conclusions based on the lack of radial segregation of the 
Blue Straggler Star (ESS) population"^. This justifies the 
simplicity of the second term on the right hand side of Equa- 
tion 1: the surface brightness profile should trace the pro- 
jected mass density. 

Having fixed the shape of the density profile, we allow 
for variation of its normalisation (and that of the total stellar 
mass Mt) by leaving M^/ Lv as a free parameter. Bellazzini 
et al. (2012) found a best-fit value of M^jLv = 1.5 ± 0.1, 
but they argue that uncertainty in stellar evolution models 
allows for a range of 1.2 < M,/Lv ^ 1-7 in the most likely 
value. For the probability density function P, we therefore 
adopt a flat distribution between 1.2 < M,/Lv < 1.7, ta- 
pered as a Gaussian of sigma 0.1 beyond these edges. 

We investigated this model parameter space with a 
Markov-Chain Monte Carlo (MCMC) approach, using the 
Parallel Tempering technique (Gregory 2005). The Paral- 
lel Tempering procedure explores the parameter space at 
different "temperature" (i.e. smoothing) levels, which en- 
sures that the solutions do not get stuck in local maxima. 
In our implementation we used 4 parallel MCMC chains, 
allowing swaps between the chains at random intervals of 
approximately 100 steps (see Gregory 2005 for a detailed de- 
scription of this method). To ensure rapid convergence, the 
proposal step sizes of all the parameters together were re- 
fined (in lock-step) every 1000 iterations, while every 10000 
steps the proposal step size of all the parameters were re- 
fined separately. We set a target acceptance ratio of 25%, 
so that on average one in four trial steps was accepted (this 
is considered to be optimal when fitting a large number of 
parameters, see Gregory 2005). 

It is well known that systems with strong radial 
anisotropy are prone to the radial-orbit instability. Typi- 
cally, stability criteria are expressed in terms of the Fridman- 
Polyachenko-Shukhman parameter ^ = 2ri.ad/?tan, where 
Trad and Ttan = + are the radial and tangential 
components of the kinetic energy tensor, respectively (Poly- 
achenko & Shukhman 1981; Fridman & Polyachenko 1984). 
Here we find it convenient to use the parameter ^haif which 



^ The argument that justifies a constant stellar M/L with radius 
is as follows: the BSS stars, which are significantly more massive 
than any evolving cluster star, are not observed to be radially seg- 
regated in NGC 2419, i.e. they have the same radial distribution 
as single stars. This should not be observed if mass segregation 
is at work, hence 2-body relaxation is ineffective in this cluster. 
Hence the stellar mix is the same at any radius, and therefore 
the stellar mass-to-light is the same at any radius. This argument 
is valid because mass segregation is expected (and observed) in 
most classical globular clusters, where 2-body encounters are a 
relevant mechanism of energy-exchange (relaxation) . In less dense 
systems such as dwarf galaxies, the BSS population (MapcUi et 
al. 2009) would not allow us to draw this conclusion, as they are 
fully non-coUisional. 



measures the same ratio, but within the half-mass radius 
(Nipoti, Ciotti, & Londrillo 2011). In paper I we found that, 
in the absence of DM, ^haif = 1.5 can be taken to separate 
stable and unstable models of NGC 2419. This condition is 
used as a prior to reject unstable solutions. 

The resulting best fit from 100000 iterations of the low- 
est temperature chain is shown with the blue line and blue 
dots in Fig. 1, and is compared to the observational data. 
The corresponding Michie model parameter values and their 
uncertainties are listed in Table 1. The solution is fully com- 
patible with the best solutions found in Paper I and Ibata 
et al. (2011b). 



2.2 Jeans equation based models 

While the Michie model gives a statistically very good de- 
scription of the cluster (as we have shown in Paper I), one 
is nevertheless left wondering how unique this solution is. 
In that earlier contribution we showed that isotropic Michie 
models (commonly referred to as King models) cannot ac- 
count for the observations, but one could attempt at this 
stage to fit alternative models such as a Plummer, a poly- 
trope or some other model. However, Nature need not follow 
these simple analytic forms. To attempt to survey the vast 
space of more general solutions, we have developed an algo- 
rithm that takes a different approach: we will assess instead 
the statistical properties of non-parametric solutions to the 
Jeans equation. 

For a spherically-symmetric system the Jeans equation 
(Binney & Tremaine 2008): 



GM{r) 2 

= — 



d In p d In cr^ 



d In r d In r 



2 1- 



2±\ 



(2) 



relates M{r) the cumulative mass inside a radius r to the 
square of the velocity dispersion profiles in the radial (cr^) 
and tangential (cr|) directions and to the logarithmic deriva- 
tives of (7^ and of the density (p). 

The algorithm we have developed explores these func- 
tions from the cluster centre out to the tidal radius. To sam- 
ple the functions in an economical manner we found it con- 
venient to assign one third of the radial bins to probe the 
inner region within two core radii with logarithmic bins. Be- 
yond two core radii the remainder of the radial points used 
linear intervals. We set the value of the central velocity dis- 
persion in the radial direction to be a model parameter. The 
radial velocity dispersion profile is then defined by the se- 
quence Aln(7^(ri), and the density profile is defined by the 
sequence A In p(ri), where Vi are the locations of the radial 
bins. To interpolate these profiles into logarithmic deriva- 
tives, we assume that the functions can be approximated 
locally (between adjacent radial bins) by parabolas in In cr^ 
(or In p) versus In r. It is worth noting that it is highly unde- 
sirable to define the profiles directly in and p, since this 
requires extremely high accuracy to avoid huge interpolation 
errors in the logarithmic derivatives. The remaining param- 
eters for the dark matter free case, are the total cluster mass 
M* and the stellar mass to light ratio Mt/Lv- 

Our model therefore has 3 parameters unrelated to the 
chosen binning ((7^(0), A/*, I^h/Lv), plus (n — 1) parame- 
ters that define the dispersion profile and n that define the 
density profile. We found it convenient to set the number 
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Figure 1. The surface brightness profile (a) and velocity dispersion profile (b) of the best-fit models. As described in the text, the 
observed data were first fit with a pure Michie model (5 free parameters) using an MCMC algorithm that explored this parameter space. 
The blue line in (a) shows the resulting excellent fit to the surface brightness profile. The model simultaneously also fits the kinematics, 
which are also shown in blue on the right-hand panel. The big hexagonal markers show a binned representation of the kinematic data, 
presented here only for visual purposes. The likelihood calculation does not use these binned values, but rather the individual kinematic 
measurements, the radial locations of which correspond to the radial locations of the small dots. The corresponding velocity dispersion 
of the model provides the ordinate value of these points. The best two-component model consisting of a Michie model to represent the 
stellar component and a double power law model to describe the dark matter is displayed in red. These values are almost coincident 
with the single (pure stellar) Michie model solution. For this composite model, a total of nine free parameters were fitted by the MCMC 
algorithm. The green and turquoise dots show the best non-parametric solutions derived with the MCMC Jeans equation solver. The 
model in green has no DM, while the turquoise model allows for a DM component. Panel (c): the surface brightness residuals between the 
data and the models (a small shift in radius has been applied to avoid the superposition of the points). Panel (d) shows the /3 anisotropy 
profile of the models. The best non-parametric model without DM (green) was derived from 1.2 X lO'^ MCMC trials, probing a vast 
parameter space of solutions; the fact its /3 profiles is very similar to that of the pure Michie model (blue) lends strong support to the 
use of the Michie model as a good representation of this globular cluster. 
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of radial bins to n = 1 + 2™, and certainly with n = 129 
the profiles are smooth and we could reconstruct a Michie 
model very accurately from artificial data. 

The algorithm projects the 3-dimensional functions 
onto the line of sight, to produce a velocity dispersion profile 
and a surface brightness profile, from which the model like- 
lihoods are measured using the observed data (Eqn. 1). One 
further assumption is made in the calculation of the likeli- 
hood: we assume that the velocity distributions are Gaus- 
sian. This assumption is not necessary for the algorithm 
to work, and we could have used the velocity moments of 
the kinematic data directly, or indeed we could have probed 
other parametrisations of the velocity distribution. However, 
it is computationally faster and avoids the errors deriving 
from the calculation of the moments from the data. Fur- 
thermore, the assumption is supported by the data: as we 
showed in Paper I, the observed distribution is consistent 
with a Gaussian function at all radii probed. The line of 
sight velocity distributions are then convolved with the un- 
certainty distribution appropriate for each data point. 

We explore this model (260 free parameters) using a 
MCMC scheme*. As before, we use the Parallel Temper- 
ing technique with 4 chains to probe progressively more 
smoothed versions of the likelihood surface. Following Gre- 
gory (2005), this is achieved by multiplying ln(L) by a factor 
where m is the chain number. The m — I chain cor- 
responds to the lowest "temperature" simulation, which is 
the one of primary interest. 

We implemented an improvement on this method, re- 
ducing the dimensionality of the high temperature chains 
by redefining their radial profiles at only n = 1 + 2™ points. 
Cubic interpolation (in In and In p versus In r) between 
these anchor points was used to fill in intermediate values 
to maintain 129 radial bins in all chains. Although we at- 
tempted a simulation including a, m — 5 chain (i.e. with 9 
radial anchor points), it was clear that the interpolation was 
too crude to be useful. 

Some regularisation is required to tell the algorithm 
that it should avoid exploring very irregular profiles. This 
was implemented by using a mild smoothing prior on the 
profiles of the logarithmic derivatives. 

An initial simulation with 2 x lO" iterations for each 
chain was run using a Metropolis-Hastings acceptance cri- 
terion. As before, all parameters were tuned in lock-step 
every 1000 iterations, and then refined independently every 
5 X 10® iterations, so as to ensure an acceptance criterion of 
25%. Finally, the 1000 best solutions from this Metropolis- 
Hastings algorithm were used as inputs for another 2 x 10^ 
iteration run using an affine-invariant ensemble sampler al- 
gorithm. Our algorithm is based on the "stretch move" al- 
gorithm proposed by Goodman & Weare (2010), which we 
expand to use the above-described multiple "temperature" 
parallel chains. The simulation uses a population of 1000 
"walkers" (which are non- coincident parameter space posi- 
tions) which evolve from one iteration to the next, and so 



The Bayesian method we are using yields probability distribu- 
tion functions for the fitting parameters, as well as information 
on the correlations between these parameters. In this context it 
is not a problem that we have more fitting parameters than data 
points. 



sample efficiently the likelihood function. The convergence 
of the method was checked by running simulations from 3 
very different starting positions. 

While the models built from the distribution function 
(described in Section 2.1) are consistent by construction, 
this Jeans-equation based approach can yield solutions that 
do not have physically possible distribution functions. This 
is a fundamental drawback of the technique, and cannot be 
avoided. However, some unphysical models are eliminated by 
ensuring that the solutions obey the Global Density Slope- 
Anisotropy Inequality (GDSAI; Ciotti & Morganti 2010a,b), 
which is the requirement that at each radius the negative 
of the logarithmic slope of the stellar density distribution 
7, = — dlnp*/d In r cannot be less than twice the local value 
of the anisotropy parameter, 7*(r) > 2/3(r-). We remind the 
reader that the anisotropy parameter is defined as: 



/3(r) 



1 - 



2a?: 



(3) 



where ar, o-^ and a^p are, respectively, the r, i9 and com- 
ponents of the velocity-dispersion tensor. As in Section 2.1, 
we also impose the requirement that ^haif < 1.5 to ensure 
stability of the resulting models. 

In Fig. 1 we display, in green, the surface brightness pro- 
file and velocity dispersion profile of the best solution to the 
Jeans equation found with this algorithm. The anisotropy 
profile of this solution is also shown (panel d); it is strik- 
ing how closely this non-parametric solution matches the 
behaviour of the best-fit Michie model from Section 2.1 
(shown in blue). That our Jeans equation solving algorithm 
should have selected this model, which is so similar to the 
previously-fitted Michie model, out of the vast space of pos- 
sible models, lends strong support to the use of the Michie 
distribution function to represent the stellar component in 
this system. The small differences with respect to the Michie 
model fit in Fig. 1 can be attributed to the differences in the 
anisotropy profile but the driving factor is the different stel- 
lar mass to light ratio {Mt/Ly ~ 1.89). 



3 MODELS WITH DARK MATTER 

Following the layout of the previous section, we consider first 
adding dark matter to the Michie model, and then freeing 
M{r) in the Jeans equation so as to also allow for dark 
matter in those solutions. 



3.1 Distribution function based models 

We first consider the possibility that the putative dark- 
matter halo can be described by the double power law den- 
sity profile (Binney & Tremaine 2008): 



Pdm(r') = 



Pdm,0 



(4) 



where 7 and 5 are the two power-law indices characterising 
the shape of the profile, is a scale radius, and pdm,o is a 
reference density. In particular, when 7=1 and 5 = 3 we 
get the Navarro, Frenk, & White (1996, NFW) model; when 
7 = and 5 = 3 we get a profile with a fiat inner density 
distribution, which is very similar to the profile proposed by 
Burkert (1995). For combinations of 7 and 5 such as those 
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mentioned above the integrated mass of the model diverges 
at large radii, so it is useful to introduce as an additional 
parameter a truncation radius rvir such that pdm = for 
r > rvir- The DM mass contained within a radius r is 



Afdm(r) = 4n pdmir'y^dr' , 



(5) 



so the total mass of the DM halo is Mdm.vir = Mdmirvir) ■ If 
the stellar density distribution is truncated at rt, the kine- 
matics of the cluster is influenced only by the mass distri- 
bution within the stellar tidal radius rt. It follows that our 
results are not sensitive to the value of rvir, provided that 
rvir > rt, which we assume throughout the paper (in other 
words, our results are not sensitive to the halo concentra- 
tion C = rvir/rs). As a consequence, we define the quantity 
Mdm.t = Afdm(rt) (the DM mass contained within rt), which 
is more useful than Afdm.vir in order to characterise our mod- 
els. In summary, our model halo is determined by specifying 
the following four parameters: 7, 5, rg and Mdm,t. 

For the DM component, for any choice of the DM pro- 
file parameters (7 and S) the adimensional potential can be 
calculated by direct integration of equation (10) as 



H/drrr(r) = 9 



Pdm,0 f r 



— I r"^pdm(r")dr" dr' 

r/r^ ^ Jo 



In the particular cases of a NFW (7=1; 5 = 3) and a cored 
(7 = 0; 5 = 3) halo the above integral admits the analytical 
solutions: 



W^dm(r) = 9 



Pdm.O ( r, 
P.,0 



2 ln(l 



(6) 



when (7 = 1,5 = 3), and 

9 Pdnr.O fr^V 77-2(1 + 77) 111(1+^) 



Wdrrr(r) 



m,0 f ^'^\ 

.,0 \rcJ 



P.,0 



f 1 + f 



(7) 



when (7 = 0, 5 = 3). 

By incorporating this dark matter potential into the 
procedure described in §2.1, one can construct spherically 
symmetric self-consistent stellar systems, truncated at a 
tidal radius rt, in equilibrium in the presence of a DM halo 
(see also Amorisco & Evans 2011). In these models the pa- 
rameters Wo, rs/vc, fa, Pdm,o/p*,o, 7 and 5 drive the shape 
of the density and velocity dispersion profiles, while the val- 
ues of p,,o and Tc determine their normalisation to physical 
units (see Appendix). 

The stellar components of the models considered here 
are self-consistent by construction, but in principle it is not 
guaranteed that there is a positive distribution function 
generating the DM components. For a wide class of two- 
component models, a necessary condition for consistency is 
that also the DM component satisfies the GDSAI (Ciotti 
& Morganti 2010b). Clearly the specific form of the halo 
anisotropy profile is not important for the purpose of the 
present investigation, because the halo influences the ob- 
servables only through its gravitational potential. However, 
we impose that in our model the halo density never increases 
with radius (dlnpdm/dlnr < at all r, i.e. 7, 5 > 0). This 
is a necessary condition for consistency only when the halo 
is isotropic or radially anisotropic, but we impose it as a 
general constraint because we think it is an astrophysically 
motivated assumption. However, we noted that it is not 



excluded that there are self-consistent tangentially biased 
models with dlnpdm/dlnr < at some r. 

Even if consistent, some models can be excluded be- 
cause they are unstable. It has been shown that DM halo 
has a stabilising effect against the radial-orbit instability 
(Nipoti, Londrillo, & Ciotti 2002; Nipoti, Ciotti, & Londrillo 
2011), so at least some models with DM halo are expected 
to be stable even if they have ^haif > 1.5. However, a full 
stability analysis of two-component models of NGC 2419 is 
beyond the purpose of the present paper, so here we conser- 
vatively assume ^haif < 1.5 as a stability criterion for all the 
models considered in the present work. With this choice we 
should effectively exclude unstable models. 

We proceed to examine the relative likelihoods of the 
two-component models built from a power-law dark matter 
halo and a Michie model stellar component. As in §2.1, we 
use a "Parallel Tempering" MCMC algorithm with 4 chains 
to explore the parameter space. The composite model has 9 
free parameters: the stellar central density, the Michie model 
core radius, the anisotropy radius, the central potential, the 
stellar mass to light ratio, the dark matter central density, 
the DM characteristic radius and the DM inner and outer 
power-law indices. For a given set of these 9 parameters, we 
construct a simulated surface brightness profile and the line 
of sight velocity distribution as a function of radius, which 
can be compared directly to the observations. 

Several (very weak) priors were imposed on the solu- 
tions. The stellar core radius must lie within a factor of 10 
of 10.88 pc; the anisotropy radius has to lie between 0.5 and 
100 core radii; the central potential parameter has to lie 
between 2 and 20; and the stellar mass to light ratio is re- 
quired to stay within the bounds of 0.8 < Mt/Lv < 2.1 (i.e. 
within 4a of the flat region 1.2 < M^/ Lv < 1-7 stated above 
in Section 2.1). Furthermore, we impose that the character- 
istic radius of the dark matter has to lie between 25 pc and 
2.5 kpc; and that the inner and outer power law exponents 
lie in the ranges < 7 < 1.5 and 2 < 5 < 4. While certain 
black hole growth models suggest that the inner power-law 
slope may be higher than 7 = 1.5 (e.g., Quinlan, Hernquist 
& Sigurdsson 1995), this is likely to affect only the very cen- 
tral region of the cluster where our data do not constrain 
the density profile. 

Figure 1 shows (in red) the most likely model found by 
the MCMC algorithm in lO'' iterations. The solution has 
Tc — 11.4pc, ra = 1.6rc, a dark matter characteristic ra- 
dius of 28 pc, and 7 = 0.06 and S = 3.7. However, we were 
surprised to find that this model requires very little dark 
matter: Mdm = 5.0 x 10*A/q inside of 10' (the approximate 
limit of the stellar structure), which is effectively insignifi- 
cant compared to the stellar mass of 8.1 X 10^ Mq. Marginal- 
ising over all other parameters, we show the posterior distri- 
bution of A/dm in Fig. 2; clearly the data appear to favour 
models without substantial dark matter. 

To confirm this curious result, we re-ran the MCMC ex- 
periments forcing certain astrophysically-motivated param- 
eter sets: cored DM models (7 = 0) and the NFW model 
(7 = 1, 5 = 3). In an almost identical manner to the more 
general solution above, the model likelihoods are found to 
decrease strongly with increasing DM fraction. In an at- 
tempt to explore whether the MCMC algorithm was able 
to escape from solutions with low dark matter, we also per- 
formed an experiment forcing a massive DM component, by 
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Figure 2. Unshaded histogram: the posterior distribution of the 
mass of the dark component inside of 10' (the approximate size 
of the stellar structure) derived by the MCMC algorithm that fits 
a Michie model inside a dark matter halo modelled as a double 
power-law. The hatched histogram shows the posterior distribu- 
tion of the luminous component of the cluster. Thus if we impose 
a Michie model on the luminous component (which as we have 
shown in Paper I, gives an excellent fit to the observations), the 
most likely solutions for the dark matter have very little mass 
relative to the stars. It would be natural to conclude that the 
cluster needs no dark component. 



setting Mdm > 5 x 10 Mq as a strong prior. But again, the 
distribution of solutions simply clustered near this low mass 
limit. 

We conclude from this analysis that a simple Michie 
model where mass traces hght is actually preferred by the 
data over more complicated models that include the possibil- 
ity of having dark matter in the form of a double power-law 
profile. The preference for a cored distribution turns out to 
be very slight - the best models with 0.9 < 7 < 1.1 are only 
a factor of 1.4 less likely than the best model. However, so- 
lutions with some DM seem better than those without DM. 
This is a manifestation of the tension between stellar and 
dynamical mass to light ratio that was noted in Bellazzini 
et al. (2012) (and also in SoUima, Bellazzini & Lee 2012 in 
other clusters) and which is discussed further below. 



3.2 Jeans equation based models 

The paucity of dark matter implied by our MCMC survey of 
the likelihood surface of Michie models is compelling, since 
the Michie models form a fairly general family of structures, 
but one is again left wondering to what extent this con- 
clusion depends on the choices of the models selected to 
examine both the cluster and the dark matter. 

For this reason we extended our analysis to a more gen- 
eral set of spherical models where both the dark matter 



Figure 3. Unshaded histogram: the posterior distribution of the 
mass of the dark component inside of 10', but this time calcu- 
lated from an MCMC algorithm that solves the Jeans equation 
in a non-parametric way. A total of 1.2 X 10^ MCMC iterations 
are presented. The difference to Figure 2 is striking: with the 
freedom from analytic form, the preferred solutions possess a non- 
negligible dark matter component. The stellar mass of the models 
(hatched histogram) is approximately half that of the dark mat- 
ter. 



density profile and the stellar density profile are free non- 
parametric profiles. We reuse the non-parametric MCMC 
tool described in §2.2 which solves the Jeans equation, find- 
ing solutions that are consistent with the data and our pri- 
ors. 

However, given that we now allow for the possibility 
of dark matter, we require some additional information to 
close equation 2. To this end we have chosen to free the 
tangential velocity dispersion profile (in addition to cr^(O), 
M, , Mt/Lv, and the profiles of A\np and Alna^)- We 
found it convenient to tabulate this profile as A In crj^, adding 
a further n — 1 parameters to the MCMC fit, plus (7|(0). The 
model thus has a total of 389 free parameters. 

The MCMC simulations were run in an identical man- 
ner to those described in §3.2, using again a hierarchy of 
4 different "temperature" chains, each with progressively 
lower spatial resolution. As before, after initial runs with 
2 X 10® iterations with a Metropolis-Hastings acceptance 
algorithm, we switched to using a population of affine- 
invariant ensemble samplers and reran the simulation for 
a further 2 x 10® iterations. 

While the models built from the distribution function 
(described in Section 3.1) are consistent by construction (at 
least for as far as the stellar component is concerned), we 
cannot exclude that some of the two-component models con- 
sidered are not generated by positive distribution functions. 
However, as in Section 2.2 we impose at least some nec- 
essary conditions for consistency. We additionally imposed 
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Figure 4. Panel (a): The cumulative mass profile of the dark matter component of the most likely model derived by the MCMC algorithm 
using the non-parametric Jeans equation method is shown with the continuous line. The dotted line shows the cumulative mass profile 
of the stellar component of the same solution. To first approximation the dark matter profile appears fairly similar to the stellar profile, 
and stands in stark contrast to cosmological expectations: the dashed line shows an NFW profile with M^ij, = IO'^Mq, concentration 
c = 23.5 and r^ir = 4.5 kpc. The dot-dashed line shows an attempt to mimic the continuous line with an NFW profile; this however 
requires an extremely high concentration (c = 250), even if r^ir = 13' = 0.33 kpc. Panel (b): the corresponding density profiles. 



the prior that the dark matter density profile decreases with 
radius and that its logarithmic slope becomes increasingly 
negative with increasing radius, i.e. if the radial locations ri 
and r2 are such that ri < r2, the logarithmic slopes at these 
locations obey the condition: 



> 



d In pdn 



dlnr 



d In pdn 



dlnr 



This property is shared by many commonly-used models 
such as the Plummer sphere, the Michie (and hence King) 
models, the double power law models if 5 > 7 (of which the 
NFW is a member), among others. 

Figure 3 shows the marginalised posterior distribution 
of the dark matter in the solutions (unshaded histogram). 
The distribution peaks at ~ 9.7 x 10^ M©, suggesting the 
presence of a significant dark matter fraction within the cen- 
tral 10'. The posterior distribution of the mass of the stellar 
component, shown with the hatched histogram, reveals that 
roughly half the amount in stars as dark matter appears to 
be present in the solutions. At first sight it would appear 
that we are facing a simple degeneracy between the total 
mass of the cluster and the mass to hght ratio of the stel- 
lar population. This suspicion is reinforced by Figure 4a, 
which shows the cumulative mass profile in the most likely 
solution: the dark matter (solid line) and stars (dotted line) 
follow each other fairly closely (but see also the density pro- 
file in panel b) . This best solution has a stellar mass to hght 
ratio Mt/Lv = 1.20, within the acceptable range found by 
Bellazzini et al. (2012). This would therefore suggest that 
the "dark matter" component is simply the result of the al- 




r [orcmin] 

Figure 5. For the most likely non-parametric solution to the 
Jeans equation (shown previously in Figure 4), we calculate the 
stellar mass to light ratio required to eliminate the need for a 
dark matter component. The region between the two dashed lines 
marks the maximum (4a-) range consistent with the analysis of 
Bellazzini et al. (2012). 
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Figure 6. The most likely Jeans equation solution of NGC 2419 using the new non-parametric MCMC technique allowing for both a 
stellar and a dark matter component. Panels (a) and (b) show the fit to the observed surface brightness and kinematics, respectively. The 
data are reproduced from Figure 1. The small green dots in (a) display the excellent model fit. The non-parametric velocity dispersion 
profiles in the radial and tangential directions {ar and a^) along with their line of sight projection (ui^os) ^'I'e displayed alongside the 
velocity dispersion estimates. We stress again that the individual star velocities were used in the likelihood calculation. The remaining 
profiles required for the Jeans equation solver are displayed in the bottom panels: the stellar density and its derivative are shown in (c), 
and the velocity dispersion derivatives in (d). For reference, we also show in panel (d) the anisotropy profile derived from the velocity 
dispersion profiles shown in (b). Note that the logarithmic derivative profiles are significantly noisier than the physical profiles. This 
model has a stellar mass to light ratio M«/Lv = 1.20, which is within the acceptable range reported by Bellazzini et al. (2012). 



gorithm attempting to force a constant stellar M* /Lv over 
the whole radial range of the cluster. 

However, this option does not appear to be possible, 
since it would require a M^/Lv profile that is inconsistent 
with observations and other astrophysical priors, as we show 
in Figure 5. Here we have calculated the stellar M^/Lv that 
would be required to account for the total density at a given 



radius. The constraints on the uniformity of the radial profile 
of the M*/Lv (Bellazzini et al. 2012), as well as the implau- 
sible Mt/Lv values beyond ~ 5' make it unlikely that the 
dark matter component in our Jeans equation solutions are 
simply due to erroneous assumptions in the M^jLv profile. 

The best Jeans equation solution is displayed in Fig. 6; 
the panels also illustrate the working of the method by 



© 2012 RAS, MNRAS 000, 000-000 



10 Ibata et al. 



showing the non-parametric profiles that are refined at each 
Markov-Chain Monte Carlo iteration. As can be seen, the 
solution gives an almost perfect reproduction of the surface 
brightness data, and interestingly, is able to fit the full range 
of the kinematic data, giving a higher (and better) central 
velocity dispersion than any of the Michie models studied 
here or in Paper I. Beyond the central ~ 1', the model is 
highly radially anisotropic^. 



4 DISCUSSION 

With current data the answer to the question of whether 
the globular cluster NGC 2419 contains dark matter or not 
depends on the model (and hence one's priors) used to anal- 
yse the structure. The choice of using a Michie model for the 
stars plus a double power law dark matter halo, is appropri- 
ate for those who value highly dynamically self-consistent 
methods and who adhere strongly to cosmological expecta- 
tions. The most likely solution of this composite model had 
Mdm{< 10') = 5 X 10*Mq, which would imply in standard 
ACDM cosmology (following Navarro & Steinmetz 2000) a 
virial mass of ~ 10^ Mq, roughly an order of magnitude 
lower than the mass of the stellar component in the cluster. 
The upper 99% confidence limit for the dark mass inside 
of 10' according to the probability distribution function in 
Fig. 2 is 7.2 x lO^M©, and this corresponds to a virial mass 
of ~ 4 X 10® Mq. Such a minuscule halo is some two orders 
of magnitude lower than expectations for the halo mass of 
dwarf spheroidal galaxies (Walker 2012), which have similar 
stellar content. Thus with this composite model the data 
prefers virtually no dark matter and one can rule out the 
possibility that the cluster is embedded at the centre of a 
dwarf-spheroidal- like mini-halo. 

The issue of whether the cluster is located at the centre 
of the putative DM halo is important for the validity of the 
analysis presented here. However, the case should in princi- 
ple be rather similar to that of the globular cluster M54 in 
the Sagittarius dSph, where the dynamical friction timescale 
for this similarly-massive globular cluster and its host halo 
to coincide is estimated to be ~ 1.5 Gyr (Bellazzini et al. 
2008). The expectation therefore is that NGC 2419 should 
reside at the DM halo centre, if it possesses one. 

The alternative analysis, using our new Jeans equation 
solver yields a slightly different answer. The strength of this 
method is that it does not require one to force a particu- 
lar analytic functional form on the data. The disadvantage 
however, is that the resulting models are simply solutions 
to the Jeans equation, and so many of them are likely to 
be unphysical (corresponding to a non-physical distribution 
function) or possibly unstable, but examining the (~ 10^) 
models individually in the detail necessary to ascertain their 
validity would be an immense task, far beyond the scope of 
this contribution. 

® Recently, Zocchi, Bertin, & Varri (2012) have studied 
NGC 2419 (among other clusters) using so-called f^"^ models, 
which can represent violently relaxed systems. In their analy- 
sis (which assumes that no dark matter is present), they also 
manage to produce a steep central velocity dispersion profile (see 
their Fig. 3). Note however that their best-fit model has 

Mf/Ly = 2.4 and is quite anisotropic (§ = 1.77). 



Nevertheless, taking the distribution in Figure 3 at face 
value gives a 99% upper bound to the dark matter inside of 
10' of 1.1 X 10® M0, which is somewhat higher than the previ- 
ous method. The strong difference is that the Jeans equation 
solution has a 99% lower limit at 7.5 x 10* M©, effectively 
demanding more dark matter than stars within 10'. The dif- 
ference with respect to the DM-free models is that that the 
solutions are better able to follow the central enhancement 
of the velocity dispersion. In order to confirm this suspicion, 
we re-ran the Jeans equation simulations from §3.2 (that in- 
clude dark matter), this time ignoring the most central 25 
radial velocity measurements. The result was a solution very 
similar to that found in §2.2, with only a small dark mat- 
ter fraction (mean total dark matter mass of 1.6 x IO'^Mq 
inside of 10'). This experiment shows that the difference in 
total dark matter content between the distribution function 
model solution of §3.1 and the Jeans equation solution of 
§3.2 is due to the latter attempting to fit better the central 
velocity dispersion, and that it is those central stars that are 
driving the requirement for a (small) dark matter fraction. 
The MCMC algorithm accommodates this mass by chang- 
ing primarily the anisotropy profile, which can be seen to be 
negative at r < 0.5' in Figure 6d (green dots). 

It is interesting to consider which model is preferable: 
the best composite Michie plus double power-law model, 
or the best model resulting from the Jeans equation analy- 
sis? The traditional frequentist approach is unfortunately 
not of much use due to the large number of parameters in 
the Jeans equation solver (we would simply find a negative 
reduced-x^, indicating that we have too many parameters). 
However, answering this question properly with Bayesian 
evidence is very challenging because the parameters used 
in the two analysis methods are very different. Indeed, it 
would require integrating the models over all parameter val- 
ues, which would be extremely computationally costly, and 
is beyond the scope of the present work. 

In this paper we have used kinematics measurements of 
NGC 2419 to test the hypothesis that this globular cluster 
is embedded in a massive DM halo. Another dark compo- 
nent that can in principle be detected via kinematics mea- 
surements is a central intermediate mass black hole (Baum- 
gardt, Makino, & Hut 2005). However, the current kinematic 
sample does not allow one to set significant constraints on 
the mass of this putative central black hole, because it af- 
fects only the very central ( ^ 0'!l) velocity-dispersion pro- 
file, which is barely probed by current observations. 

4.1 Comparison with previous work 

Baumgardt et al. (2009) assume that the putative halo 
around NGC 2419 has a NFW profile (with = 500 pc) 
and use their data to limit the dark mass inside of 500 pc 
to Mduiir < 500 pc) < IO^Mq. Given their model param- 
eters, this corresponds to Mi^-,{r < 10') < 3.9 x 10® M©. 
In the work of Conroy, Loeb, & Spergel (2011), they also 
assume that the halo has a NFW profile (with rs — 250 pc 
and Tvir = 1 kpc) and conclude that M^mir < Ikpc) < 
10® Mq, which, for their model parameters, gives a very tight 
Afdm(r < 10') < 2.4 X 10'*Mq. Our results for the two com- 
ponent models are in reasonable agreement with these pre- 
vious results. The limits provided above, are more stringent 
than those of Baumgardt et al. (2009), and unlike Conroy, 
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Loeb, & Spergel (2011), our method does not require the 
cluster to have suffered strong two-body relaxation, which as 
we have mentioned, should not have occurred in NGC 2419. 



5 CONCLUSIONS 

We have undertaken an analysis of current data on the kine- 
matics, structure and stellar mass to light ratio of the Galac- 
tic globular cluster NGC 2419. This distant halo cluster is a 
very interesting specimen to study, since it is massive and as 
such may have been able to maintain a residual dark mat- 
ter halo over the course of its assimilation into the Milky 
Way. Thus NGC 2419 allows us to examine the possibility 
of whether (some) globular clusters formed along with (their 
now disrupted?) dwarf galaxy hosts within mini dark matter 
halos. 

Our earlier analysis (Paper I) showed that a simple 
Michie model gives an excellent fit to the state-of-the-art 
measurements of the kinematics and surface brightness pro- 
file of the cluster presented in that contribution. Here we 
extended that study to investigate whether composite mod- 
els with an additional dark matter component in the form 
of a generic double power-law could also be accommodated. 
To this end we implemented an MCMC parameter search 
for solutions of this composite system. The Michie distribu- 
tion function is found to be highly restrictive; while it is as 
yet not possible to completely rule out the presence of dark 
matter in NGC 2419, the preferred dark matter fraction is 
6% of the stellar mass within 10' (the approximate limit of 
the cluster), and there is virtually no room for astrophysi- 
cally interesting amounts of dark matter that could hint at 
the formation of clusters within dark matter halos. 

In a major effort to assess the model-dependence of 
these results we developed a new Markov-Chain Monte Carlo 
algorithm that solves the spherical Jeans equation using 
non-parametric profiles for the density and velocity disper- 
sion profiles. A very high number of free parameters is re- 
quired to follow the gradients with sufficient resolution. The 
code uses a hierarchy of resolutions so as to probe the pa- 
rameter space as efficiently as possible. These ideas (and the 
software) should be applicable in the future to analyse other 
approximately spherical systems of interest, such as dwarf 
galaxies, elliptical galaxies, or galaxy clusters. 

The more general Jeans equation solutions prefer a 
small dark matter fraction, approximately twice the stel- 
lar mass content, and at face value rule out the possibility 
that there is no dark matter. This last inference depends, 
however, on the validity of the conclusion by Bellazzini et 
al. (2012) that the stellar mass to light ratio is invariant with 
radius. For the most likely solution, the inferred dark matter 
profile has a very different shape to an NFW model, so it is 
possible that this putative dark component, if real, is related 
to unexpected stellar remnants (the expected remnants were 
taken into account by Bellazzini et al. 2012) rather than cos- 
mological dark matter (note however, that the resolution of 
present-day simulations that include the relevant baryonic 
physics do not allow reliable predictions for the profile of 
dark matter mini-halos of this mass) . 

In their analysis, Bellazzini et al. (2012) found that as- 
suming a 100% retention of all the dark remnants would 
have changed the best fit M,/Lv from 1.49 to 1.6. Hence in 



a standard scenario, changing the fraction of remnant stars 
will not have a sufficiently large effect to account for Fig. 5. 
On the other hand, as is pointed out in that contribution, 
all the currently available (and highly uncertain) models in- 
voked to explain the multiple populations in GCs envisage 
(a) a progenitor of the present-day cluster with a baryonic 
mass a factor ^ 10 larger than today, (b) a preferential loss 
of first generation stars with respect to the second genera- 
tion, (c) a strong difference in the radial distribution of stars 
of the two generations (the second one being more concen- 
trated, as is generally observed; see Lardo et al. 2011). This 
complex star formation history and strong variations in the 
cluster potential must have occurred in the first < 500 Myr 
of the cluster life. It can be conceived that at this stage both 
(i) an anomalous production and retention of dark remnants 
occurred, and (ii) a radial distribution of dark remnants dif- 
ferent from that of the light was produced and then remained 
un-modified in that non-coUisional environment. 

Therefore if the observed dark matter is of baryonic 
origin it is likely that it has something to do with the origin 
and early enrichment of the multiple populations (Cohen, 
Huang & Kirby 2011; Mucciarelli et al. 2012) inhabiting 
NGC2419. On the other hand one may imagine that we see 
the small remnant of a larger non-baryonic cosmological halo 
that was nearly completely stripped away from the cluster at 
early times (see, e.g. Pefiarrubia, Navarro & McConnachie 
2008). During the first approach to its peri-galacticon, it 
may have provided the potential well that allowed the cluster 
to retain ejecta from first generation stars. 

Thus all the models we have explored turn out to have 
very little dark matter in their central regions, and we rule 
out any connection with mini dark matter halos of the mass 
that is associated with those in which dwarf galaxies are 
embedded. This conclusion, though consistent with earlier 
work on this particularly promising cluster, and consistent 
also with a large literature on more nearby systems (but 
which on the whole are clearly not simple isolated systems), 
remains puzzling. The connection between globular clusters 
and dwarf galaxy remnants is becoming more and more com- 
pelling (Mackey et al. 2010), so it is natural to expect that 
some clusters retained some of the dark matter of their for- 
mer host. It is possible that this problem will only be re- 
solved conclusively by obtaining a larger sample of distant 
halo clusters, but that will require observing further afield 
in M31, and to be feasible will require the next generation 
of instruments. 
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Table 1. Model parameters and derived quantities for the four families of models investigated in this contribution. Where possible, we 
report the mean value and root mean square dispersion of the MCMC solutions. Due to the construetion of our Jeans equation algorithm, 
the stellar eentral density is expensive to calculate at every MCMC iteration, so we report only the values of the most likely solution. 
Furthermore, the posterior probability distribution for the dark mass in Fig. 2 is highly skewed, so we only report the most likely value 
for that model. 



Michie Jeans without DM Michie + DM Jeans with DM 



stellar central density {Mq pc'^) 57.3 ± 4.8 74.4 (best value) 59.2 ± 10 66.9 (best value) 

stellar core radius ( pc) 12.6 ± 0.6 ... 10.8 ± 0.5 

anisotropy radius ( pc) 14.2 ± 2.3 ... 18.3 ± 2.7 

central potential 3.6 ± 0.4 ... 4.7±1.1 

stellar M/L 1.72 ±0.11 1.89 ±0.07 1.62 ±0.17 1.15 ±0.12 

Stellar mass (IO^Mq) 8.1 ±0.2 9.0 ±0.1 7.6 ± 0.8 5.5 ±0.6 

Dark mass (lO^M©) 0.5 (best value) 9.7 ±0.8 



APENDIX A 

The distribution function proposed by Michie (1963) has the 
following phase-space form: 



f{E,L) = fo exp - 



exp 



f{r,Vy,vt) = /o exp 



exp 
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2 I 2 
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(8) 



where E and L axe the energy and angular momentum per 
unit mass, respectively, ip is the effective potential defined 
as the difference between the cluster potential (including the 
contribution of both stars and the dark matter — the latter 
component will be added in Section 3.1) at a given radius r 
and the potential at the cluster tidal radius {ip = 4' ~ 4't), 
fo is a normalization coefHcient factor, ax is a reference 
velocity dispersion, Vt and are the radial and tangential 
component of the velocity and ra is the anisotropy radius: 
the velocity distribution tends to be isotropic at r < ra and 
radially anisotropic at r > ra. 

The density and the radial and tangential components 
of the velocity dispersion can be derived by integrating the 
above distribution function: 
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It is convenient to work with the dimensionless quantities 
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where p«,o is the central cluster density and 
r. = (^^^ 

is the core radius (King 1966). The above equations can be 
then written as 
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which can be solved once the potential W{r) is known. This 
last quantity can be written as the sum of the stellar and 
DM contributions W = Ty* + W^dm- Both contributions must 
obey independently the Poisson equation 
72 



d^W, 2dm 
df2 f dr 



'9p, . 



(10) 



Here, the subscript i refers to either stars or DM. For the 
stellar component we leave Wo as a free parameter and solve 
equation (10) adopting the boundary conditions at the cen- 
ter 

W.,0 =Wo- Wim,0 , 

dW,,o 



dr 



■0. 



Equations (9) and (10) can be integrated numerically to 
derive the three-dimensional density profile, and the radial 
and tangential velocity dispersion profiles. Finally, the above 
profiles have been projected on the plane of the sky to obtain 
the surface mass density 



E.(i?) 



p,rdf 
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(11) 
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and the line-of-sight velocity dispersion 



i?2 



(12) 



As described in Section 2.1, the full likelihood analysis 
performed here requires for each model the knowledge of the 
distribution of projected velocities (not only its second-order 
moment). This function can be calculated as 



fpiv,R) = 



p* (r)r 



(13) 



— 2il> — v'^—v'^ 

X / f{r,v'^,v't) dwydw^dr , 

Jo 

adopting the following change of variables 
v'y. — V sin a + Vx cos a, 

III • \2 , 2il/2 

Vt — [(v cosa — Vx sma) +Vy\' , 
where a = arccos(i?/r). 
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